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ABSTRACT 


An investigation of inertial coupling and its contribution to wing rock in the 
F-14A aircraft has been conducted. Wind tunnel data was used to obtain the 
stability parameters for angles of attack from zero to 25 degrees, after which 
linear and nonlinear analyses of the equations of motion were completed. The 
linearized analysis of the uncoupled longitudinal and lateral-directional equations 
was included to provide a baseline for comparison with the fully coupled, 
nonlinear equations. In both cases, the equations of motion were solved 
numerically and time history traces produced to illustrate aircraft response. 
Results indicate that a stable short period mode can feed damping energy into an 
unstable dutch roll mode via the coupling of the equations to produce a stable 
limit cycle very similar to those experienced in the aircraft. Numerous 


suggestions for follow on research are presented. 
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I, INTRODUCTION 


Much has been written on the subject of linearized aircraft flight mechanics. 
Under the assumption of small perturbations, the equations of motion 
“representing rigid body aircraft dynamics can be reduced to two sets of 
uncoupled, linear differential equations with constant coefficients which can be 
Teadily solved to yield characteristic frequencies, damping ratios, and 
- amplitude/phase relationships of the aircraft's "natural modes", as well as time . 
history traces of aircraft response to these modes due either to initial conditions 
or contro] inputs. The approximations made in th. analysis of the linearized 
equations are quite good as long as the aircraft motion is analyzed within the 
boundaries of linear behavior (i.e., small angles). Under the assumptions made 
in this relatively benign region, the longitudinal aircraft motions are isolated 
from the lateral-directional motions. Analytical results from linearized theory 
agree well with flight tests in this regime, with negligible coupling occurring in 
flight. As good as these approximations are, they break down completely at high 
angle of attack, sideslip angle or high angular rates and fail to provide accurate 
results when analyzed in this manner. Thus, the full set of nonlinear, coupled 
equations of motion must be analyzed to obtain adequate results. Whereas the 
motion and inertially related nonlinear terms were neglected in small 
perturbation theory, these terms provide coupling between the longitudinal and 
lateral-directiona! dynamics in the nonlinear analysis, altering the overall 
dynamic response. 

The differences between linear and nonl‘near flight dynamics described 


above form the foundation for this study. In particular, this analysis probes into 














the nonlinear phenomena known as wing rock. Numerous studies conducted 
over recent years in an attempt to identify the specific cause(s) of wing rock have 
highlighted aerodynamic hysteresis (Ref. 1], complex asymmetric vortex 
interactions for slender deltas (Ref. 2, 3) and inertial cross coupling (Ref. 4] as 
potential candidates for the wing rock motion. The validity of these studies is 
certainly not under question; however, not one of them can be chosen as the 
specific cause for wing rock in any aircraft without careful consideration due to 
the highly configuration-dependent nature of the problem. It is with this in mind 
that this study investigated the effects of inertial cross coupling as a contributor 
to the wing rock motion in the F-14A aircraft. Although wing rock has been 
reported in many other aircraft (e.g., A-4, F-4, F-5, (Ref. 1]; HP-115 [Ref. 2); 
T-38 (Ref. 4); F-15 (Ref. 5]; F/A-18, X-29 (Ref. 6]) this aircraft was modeled for 
a variety of reasons. Perhaps most important among these reasons was the 
rapidly growing emphasis on high angle of attack research and technology, along 
with its potential impact on the future combat capability of the F-14. The very 
existence of aircraft currently displaying excellent high angle of attack 
performance such as the F-16, F/A-18, and Mig-29 merely suggest that future 
designs will display even greater capabiiity in this regime. As the F-14 enters its 
third decade of service, it carries with it a battle-proven record of superior 
performance. But the realities of reduced military spending and a huge federal 
deficit make it all too clear that every weapon system in the inventory will be 
used to the very end of its service life (and as we are witnessing with A-6's, at 
times beyond the projected service life by sending aircraft through rework 
facilities for structural repair and/or enhancement; a less time consuming and 


less expensive alternative to development, test and acceptance of new designs). 

















Thus, if trends in the defense contracting industry and weapons procurement 


~ ‘branches of the armed services continue to proceed as they have in the 80's, it is 


not unreasonable to assume that the F-14 will remain in service for the 
foreseeable future, meaning at least into the 21S¢ century. The F-14 will 
undoubtedly face aircraft with superior high angle of attack capability and may 
be forced, depending on rules of engagement, current tactics, etc., to fight these 
aircraft in the high AOA regime. Even with the current upgrades to engines and 
avionics, the basic airframe of the F-14 will remain the same. Therefore, the 
dynamic behavior of the aircraft from a stability standpoint will also remain 
unchanged unless alterations are made to the aircraft's stability augmentation 
system. If an F-14 experiences wing rock while engaged in an aerial encounter, 
especially at very high angle of attack, the only current recourse is to neutralize 
lateral and directional controls and momentarily reduce the angle of attack. 
Obviously, this course of action may compromise the aircraft and prove to be 
unrecoverable for the crew. If it is within our means to understand this 
phenomenon, then we are that much closer to providing and implementing a 
solution in the form of control laws to avoid it. 

Another important reason for studying the F-14 was the relatively simple 
flight control system incorporated into its design. As opposed to many of the 
newest high performance aircraft, the F-14 flight control system uses only pilot 
input via mechanical linkages and a three axis stability augmentation system 
(SAS) to control the aircraft during normal maneuvering. The aircraft can be 
flown freely within its envelope without pitch SAS and is restricted to roll SAS 
off above 15 units AOA for subsonic maneuvering. Thus, only the function of 


the yaw SAS, the most critical of the three, was of concern to this study. 
































Although one design feature of the yaw SAS was to provide increased directional 
‘Stability at high angle of attack, the aircraft can still be maneuvered normally 
without it "if extra care is taken to control yaw excursions with rudder." (Ref. 
7:p. IV-11-1] Therefore, conducting an analysis of the stick fixed stability 
characteristics of this aircraft without knowledge of the details of its stability 
augmentation system is not unreasonable. For aircraft with "fly-by-wire" flight 
control systems, the results of a similar stability analysis may vary by a wider 
margin (from actual aircraft response) due to the influence of the flight control 
computer on the control surfaces in response to changing flight conditions, even 


with no pilot input. 


The F-14 also proved to be an excellent choice due to the availability of an 


extensive aerodynamic data base. Lastly, an abundance of pilot reports from 
fellow aviators at NPS was readily available. 

In an effort to ultimately obtain numerically derived time history traces 
illustrating wing rock for the F-14, the sections that follow provide a brief 
background on F-14 high AOA flight characteristics, the aerodynamic data base 
from which the aircraft stability derivatives were obtained and a summary of the 
equations of motion. The study continues by describing the numerical 
procedures used, then highlights the linearized and nonlinearized results. A 
number of recommendations for related research are presented following the 


conclusions. 
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Il, F-14A FLIGHT CHARACTERISTICS REVIEW 


A. NATOPS FLIGHT MANUAL REVIEW 

The F-14A NATOPS Flight Manval (Ref. 7] mentions wing rock and 
coupling tendencies in the Flight Characteristics section of the manual, The 
following quotations from the manual are provided to familiarize the reader with 


the nonlinear dynamics of this aircraft. 


Coupling occurs when motions in more than one axis interact. The F- 
14A, like all high-performance airplanes capable of producing high rated, 
multiple axes motion, is susceptible to coupling. High rate, multiple axes 
motions particularly at high AOA can produce violent coupled departures. 
(Ref. 7:p. 1V-11-18] 


Although the maneuver slats increase the severity of the wing rock 
between 20 and 28 units AOA, overall departure resistance of the aircraft is 
greatly improved. The wing rock may be damped with rudders, but greater 
difficulty may be encountered with maneuver slats, especially at low 
airspeeds. If this occurs, the wing rock may be damped by neutralizing the 
lateral and the directional controls and momentarily reducing AOA to 
below 20 units. (Ref. 7:p. IV-11-7] 


At 20 to 28 units AOA, reduced directional stability is apparent, and 
even small control inputs will cause yaw oscillations that, if unchecked, can 
produce a mild wing rock (+/- 10 to 15 degre: s).... Wing rock at 20 to 28 
units AOA will be more severe (+/- 25 degrees) and more difficult to damp 
with maneuver slats extended (due to increcsed dihedral effect). [Ref. 7:p. 
IV-11-11) 








In the takeoff and landing configuration, wing rock is also experienced at 
high AOA during approach to stall, as indicated by the following excerpt from 
the NATOPS. 


At 25 units AOA divergent wing rock and yaw excursions define the 
stall. Sideslip angle may reach 25 degrees, and bank angle 90 degrees 
within 6 seconds if AOA is not lowered.... Extending. the speed brakes ... 
improves directional stability significantly, reducing the wing rock and yaw 
tendency at 25 units AOA.. Stall approaches should not be continued beyond 

. the first indication of wing rock. When wing rock occurs, the nose should 
be lowered and no attempt should be made to counter the wing rock with 
lateral stick or rudder. [Ref. 7:p. IV-11-18] 


B. FLIGHT TEST TIME HISTORIES 

To further illustrate the characteristics of wing rock for this aircraft, several 
time history traces from two actual flight tests are shown in Figures 1a-1d and 
2a-2d [Ref. 8). In both sets, it can be noted that the lateral acceleration of the 
c.g. shown in Figures 1d and 2d can be converted into a sideslip perturbation. 

The first set of traces, Figures la-1d, was obtained from a flight test 
conducted at approximately 20,000 ft., Mach .32, with maneuver flaps extended, 
wings swept fully forward and the SAS off. The initial angle of attack (ARI 
ALPHA) from the plot was 14 degrees. This angle corresponded to 


approximately 12 degrees true AOA via the conversion 


AOArue = .8122*AOA aR] + .7971 (01) 





fA ee ae a 


(Ref. 8] and 17 units indicated AOA via the conversion 
AOAynits = 1.0989 *(AOAtrue + 3.01) (02) 
which is valid for the F-14 operating at Mach number less than .4 [Ref. 9]. 


As the time history began, the aircraft was decelerating and most likely 


experiencing light buffet (note the rapid variations in normal acceleration, 


‘Figure 1a). In order to maintain altitude, the pilot commanded aft stick. This 


input resulted in corresponding increases in stabilator position, angle of attack ~ 


and pitch attitude (Figure 1b). Around 15 seconds and at the equivalent of 
approximately 18.5 units AOA, very slight roll rate oscillations began which 
resulted in a 20 degree left wing down attitude after about 28 seconds and 
coupled into small angle of attack and sideslip perturbations. The pilot 
responded with a slight right stick correction towards wings level, at which time 
larger roll rate and roll amplitude oscillations began to develop, at around 30-35 
seconds (Figure lc). The angle of attack at the point where pronounced wing 
rock oscillations first appeared was the equivalent of 20.5 units. At the same 
time, larger sideslip and yaw rate oscillations began to appear (Figure 1d), as 
well as a rapid nose down pitch rate. As the left wing down oscillations were 
arrested and the mean value of the wing oscillation retumed towards wings level, 
the lateral stick input was taken out. With the decrease in pitch attitude came a 
loss in altitude and an increase in velocity. Recovery began by applying forward 


stick to reduce angle of attack at about 45 seconds, although the oscillations 


continued beyond the end of the trace at 60 seconds. 
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Figure Ic. Flight Test Time History Traces 
F-14A 1X, BUNO 157991, FLT 264, GW 55386 LB., C.G. 6.8 % MAC 
SWEEP=19 DEG, FLAPS=11 DEG, ALT=20200 FT, MACH .32, SAS OFF 
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Figure Id. Flight Test Time History Traces 
F-14A 1X, BUNO 157991, Fi.T 264, GW 55386 LB., C.G. 6.8 % MAC 
SWEEP=19 DEG, FLAPS=11 DEG, ALT=20200 FT, MACH .32, SAS OFF 
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The second set of traces, Figures 2a-2d, was taken from thc same test flight 
at a nearly identical flight condition, except that the altitude was approximately 
19,000 ft. The initial angle of attack (ARI ALPHA) from the plot was 17 
degrees. This corresponded to approximately 14.5 degrees true AOA and 19.5 
units as seen in the cockpit. In this demonstration, the aircraft began by quickly 
decelerating and increasing AOA to 19 degrees ARI (approximately 16.5 degrees 
true and 21 units, Figures 2a and 2b). The immediate presence of roll and roll 
rate oscillations without lateral control inputs is evident in Figure 2c, as well as 
sideslip and yaw rate oscillations in Figure 2d. The roll amplitude eventually 
reached approximately +/- 45 degrees with a period of about 5 seconds. It is 
evident that frequency doubling in the AOA perturbation occurred, as the period 
was about 2-2.5 seconds. The AOA reached approximately 25.5 units at its 
highest point. Also, a pronounced decrease in pitch angle was evident once small 
AOA perturbations appeared at around 10 seconds. Furthermore, loss of altitude 
appeared as the roll angle oscillations built beyond +/- 30 degrees from the mean 
value at around 25 seconds. 

All of these illustrations clearly indicate that the motions which occur at high 
angle of attack are quite complex and that nonlinear coupling occurs between the 


longitudinal and lateral-directional equations of motion. 
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Figure 2b. Flight Test Time History Traces 
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Figure 2c, Flight Test Time History Traces 
F-14A 1X, BUNO 157991, FLT 264, GW 55345 LB., C.G. 7.1 % MAC 
SWEEP=19 DEG, FLAPS=11 DEG, ALT=1898£ FY, MACH .32, SAS OFF 
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Figure 2d. Flight Test Time History Traces 
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Hil. A DESCRIPTION OF THE AERODYNAMIC DATA BASE 


The stability derivatives needed to obtain the coefficients in the equations of 
motion for the F-14 came from an aerodynamic data base provided by the 
- Grumman Corporation, Bethpage, New York [Ref. 9]. The data base consisted 
of aerodynamic and stability coefficients from wind tunnel and spin tunnel tests 
conducted on a flight test aircraft at the NASA Langley Research Facility. The 
tabulated data represented two different flight regimes; a low speed regime 
where Mach number remained at or below M=.6 and a high speed regime for 
Mach number greater than M=.6. Typically, the low speed data appeared in the 
tables as a function of angle of attack and sideslip, although some stability 
derivative coefficients were dependent solely on AOA, while others were 
dependent upon AOA, sideslip and control surface deflection. All data in the 
tables was referenced to the body axis coordinate system (x - axis corresponding 
to the aircraft's longitudinal axis). Therefore, coordinate transformations 
through the angle of attack were required to obtain coefficients related to lift and 
drag. Additionally, crossplots, curve fits and a number of other techniques were 
employed to obtain all of the necessary information for the stability analysis. A 
detailed description of the techniques used to obtain the stability parameters from 
the wind tunnel data appears in Appendix A. 

It is very important to remember the distinction between true angle of attack 
expressed in degrees and in units. The conversion between degrees true and 
units for the F-14 at Mach number less than .4 was given in equation (02). The 


net effect of this conversion is to subtract five from the indicated AOA expressed 


in units to get AOA in degrees (i.e., 25 units = 20 degrees). This distinction 








should be kept in mind when comparing the NATOPS flight characteristic s to the 
results of this study, which indicate aircraft response at various trugé angles of 


~ attack expressed in degrees. 
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IV. EQUATIONS OF MOTION DEVELOPMENT 


The equations of motion describing the dynamic behavior of an aircraft are 
typically developed in a rotating reference frame fixed to the aircraft's body 
. .axes, The aircraft is assumed to be a rigid body, such that aeroclastic and : ie 
gytoscopic effects are neglected. Additionally, the variation in aircraft velocity 
is assumed to be negligible for the purposes of this study. This approximation 
allows for the removal of the X force equation and the perturbation velocity 
from the longitudinal equations of motion. After the assumption of small - 
perturbations is made and the linearization of the equations is completed, they 


can be conveniently expressed in state space format as shown below. 
Xiong = Along Xlong (03) 
and 
Xlat-dir = Alat-dir Xlat-dir (04) 


where 


Xlong = [0 9 8 } (05) 








U 
: Along = Mo! Mg 0 (06) 
0 
1 
Ma’ = Ma +75 (Mg Za) (07) 
Mg' = Mg+Mg (08) 
: T 
Xlat-dir = Bpo | (09) 
YB gcos8Q 
yo gg 
L 0 L 
Meat? / (10) 


NB Np O- Ny 


This method of describing the governing equations is compact and readily 
shows that the longitudinal and lateral-directional equations for the linearized 
case are uncoupled. It also lends itself to matrix operations for determining the 


stability characteristics of the natural modes for both longitudinal and lateral- 
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directional motion. A general description of the natural modes normally 

associated with small perturbation theory is presented in Etkin [Ref. 10]. The 

“ dimensional derivatives which make up the elements of Ajong and Ajat-dir are in 

accordance with NASA convention and are defined in McRuer, Ashkenas and 
Graham (Ref. 11). 

_ When the underlying assumptions of small perturbation theory are extended 
in order to introduce non-linear terms in the equations of motion, the linearized 

- equations shown above must be modified to account not only for the extra non- 
linear terms, but also to drop out the terms resulting from the linearization. 
Additionally, Euler angle relations are introduced so that the roll angle o 
becomes the Euler angle ® and the pitch angle © becomes the Euler angle @. 
The additional terms to be added to transform the linearized equations into the 


fully coupled, non-linear equations are shown below. 





(NL)p = po + £ (CaSp - Cao6) (11) 
ly ° Iz > 
(NL)p =" at (12) 
(NL)Q = (qSo + Co) Tg (13) 
ily 
(NL)p=—— ap (14) 
Zz 
(NL)qx = -pB + © (CgCo - Cay + Seq) (15) 
Iz is Ix 
(NL)q Ske Tp + My (NL)o (16) 
(NL)g = (Co - 1) - 89 (17) 


zi 








and the new state vector for the non-linear equation set is: 
| T 
XnL=[B p © ra q @] (18) 


It will be noted that the aerodynamic terms still retain their linear form 
while the nonlinear aspects are introduced via inertial coupling. Retaining the 
simplified forms for the aerodynamic terms is deliberate at this stage of the 
analysis in order to illustrate the influence of inertial coupling upon the ensuing 


motion response. 
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V. COMPUTATIONAL PROCEDURES 


A. LINEAR SYSTEMS OF EQU:TIONS 

_ An understanding of the equations of motion and the significance of their 
“nonlinearity is essential prior to attempting to determine the aircraft's response 
to an initial perturbation. Had the equations remained purely linear, a relatively 
simple solution would have been available by calculating the eigenvalues of the 
characteristic polynomial and each associated eigenvector. The general solution’ 
of a linear differential equation is a linear combination of all linearly 
independent solutions. Therefore, the general solution can be expressed 
explicitly as a function of time as follows (two degree of freedom system shown 


for clarity): 
(x()) =cy (x1) el £09 (x9) aoe (19) 


where {x(t)) is a column vector whose individual components represent the time 
response of each degree of freedom, A; and Az are the eigenvalues, {x1} and 
{x2} are the associated eigenvectors and cj and c2 are constants representing the 
modal participation of each of the modes in the response. In general, the 
eigenvalues are complex in nature and represent a damped, oscillatory solution. 
Utilizing the fundamentals of complex arithmetic as applied to a response with 
real physical terms, the expression may be rewritten as a combination of sines 
and cosines to remove the complex terms, thereby leaving the equation as a 


trigonometric function of time alone. Once the time domain solution is found in 
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this manner, time history traces may be generated using a short computer 
algorithm. Such an algorithm might iterate time, calculate the value of each 
component of ({x(t)} at that time, store the values in a data file and repeat the 
procedure. Plotting the data would provide the time history traces. Although 
the procedure outlined here was intended to iilustrate the ease with which a 
purely linear system can be solved, the technique described in the next section 
was used to solve both the linear and nonlinear systems of equations so that only 
one program had to be used regardless of the desired form of the solution. Even 
so, the calculation of eigenvalues, eigenvectors, natural frequencies and damping 
ratios remained important and was used extensively in the linear analysis because 
the character of the dynamic response could be determined and visualized solely 
by these parameters. The actual time history response provided an added means 


of visualizing the influence of these parameters. 


B. NONLINEAR SYSTEMS OF EQUATIONS 

A computer program was used to numerically integrate both the linear and 
nonlinear equations of motion to obtain a time history response. As mentioned 
earlier, the algorithm was capable of numerically integrating linear differential 
equations as well as nonlinear equations. Therefore, the same program was used 
to obtain either the linear or the nonlinear response, depending on the solution 
desired. 

The heart of the program was a numerical integration technique based on 
Richardson's extrapolation; an explicit, two-step numerical procedure which can 


be readily applied to first order differential equations. It is a variation of the 
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second order Runge-Kutta method. The basic concept behind the method is 
presented by Ferziger [Ref. 12] and demonstrated by the following example. 
Suppose the first derivative of each variable of interest is expressed 
explicitly as a function of known initial conditions at time t=0. By simple 
substitution, the numerical value of each derivative can be calculated and used to 
approximate the value of the variables of interest a short time later. 


Mathematically, 
(1) 
Y(n+1)=Y¥n*thyn (20) 


where h is ‘he time increment, n represents a discretized time step and the 
superscript (1) denotes that the value of y(n+1) calculated here is just the first 
estimate of the final value at time t(n+1) = tn + h. The calculation is then 
repeated using two steps, each at half of the original time increment as shown 
below. 


(2) h: 
oy ae ae is (21) 


(2) 
Y(n +1) = Y(n Mt z y(n +7) (22) 


The superscript (2) indicates that the calculated value is the second estimate 


of the final value. The estimates are then combined using Richardson's 
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_ extrapolation to obtain the final expression for the variable of interest at the new 


time. 


(2) (1) 
Y(n+1)=2 Y(n+1)* ¥(n+1) (23) 


This procedure is repeated for all of the variables at each time step to obtain 
updated values for each variable at the new time step. The value of each variable 
at the new time step is then stored in a data file for subsequent plotting. The 
time is then incremented and the first derivatives recalculated using the updated 
information from the previous time step. The process continues for a length of 
time as specified by the user, at which time the computer program terminates. 
After an investigation to determine the time step required to retain sufficient 


accuracy was conducted, the time step was set at .05 seconds. 


C. INPUTS TO COMPUTER PROGRAM 

The program mentioned in the previous section required specific inputs 
which served to identify the geometry, configuration, inertial properties and 
aerodynamic characteristics of the aircraft in its trimmed condition. While most 
of this information remained constant regardless of the angle of attack studied, 
the aerodynamic characteristics which determine the dynamic response of the 
aircraft were highly dependent upon, among other things, angle of attack. 
Therefore, to describe the character of the aircraft in the program at each 
different angle of attack, the longitudinal and lateral-directional “plants” 
corresponding to the linear response at the desired trim AOA were used. These 


“plants” are the square matrices labeled Along and Alat-dir in equations (06) and 
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~ (10) which contain the coefficients in the equations of motion. These coefficients 
are made up of specific combinations of the stability derivatives which are 
“normally determined by wind tunnel testing and post flight parameter 
identification techniques. The nonlinear terms were programmed in and selected 

as a program option if a nonlinear solution was chosen by the user. 

In addition to these items, initial conditions were specified to obtain a non- 
zero response. The choice of initial conditions was important in obtaining a 
coupled response consistent with the documented response of the aircraft. The 
choice of initial conditions will be discussed later. A complete listing of the 


computer program appears in Appendix B. 
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VI. ANALYSIS 


The ultimate goal of this analysis was to obtain multiple time history traces 
which characterized the aircraft wing rock motion, one trace corresponding to 
each of the seven pertinent degrees of freedom. In building up to this goal, a 
“linearized analysis (and, therefore, uncoupled as well) was used first to ensure 
that the aircraft's response to known stabilizing or destabilizing conditions would 
produce convergent or divergent behavior, respectively. Furthermore, the 
linear results could be used as a benchmark to compare to the nonlinear, coupled 
response once that response was determined. The study included aircraft 
response at angles of attack ranging from zero to 25 degrees so that the 
variations in the response due to the initial trim condition could be evaluated. A 
summary of the trim conditions for flight with gear up and flaps at maneuver at 


500 ft. appears below in Table 1. 


TABLE 1 
SUMMARY OF AIRCRAFT TRIM CONDITIONS 


Trim Velocity] Trim oe Mach 
iseale units ft/sec kts number 
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The geometric and inertial properties of the aircraft described in the data 
base were as follows: | 
oe SN "WT = 52000 Ib. 
C.G. = 16.2 % MAC 
Ixx = 51509 slug ft2 | nos 
‘yy = 232773 slug ft2 - | —— 
Izz = 275627 slug ft2 


Ixz= 2654 slug ft2 


It is critically important to recognize that although an aircraft is never flown 
stick fixed, a stability analysis conducted in this manner provides very useful 
information on the tendencies of the aircraft to move about its trim position once 
disturbed. The flying qualities of an aircraft are very much dependent on these 
stability characteristics. In fact, handling qualities ratings are based on the 


dynamic characteristics of the linearized modes for many aircraft. 


A. LINEAR ANALYSIS 

To begin the linear analysis, the characteristic polynomial for both the 
longitudinal and lateral-directional "plants" were solved for the eigenvalues of 
each system. The natural frequencies, modal damping and associated 
eigenvectors were also calculated. As previously mentioned, these parameters 
serve to identify the character of the dynamic response for each of the natural 
modes. A summary of the linearized systems short period and dutch roll 


characteristics is provide below in Table 2. 
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TABLE 2 
SUMMARY OF LINEARIZED DYNAMIC CHARACTERISTICS 


Short Period 


a Natural Frequency | Damping Ratio Damping Ratio 
_. .L (degrees rad/sec rad/sec 





A quick look at this table reveals a great many things about the character of 
the aircraft's dynamic response at the angles of attack studied. The following 
paragraphs provide an in depth discussion of the linearized results for the F-14A. 

1. Short Period Mode 

a. Time History Response 
The F-14 displays a stable, heavily damped short period mode at all 
of the angles of attack studied in this analysis. When excited, this mode tends to 
produce a rapidly convergent solution back to the initial trim condition. The 
short period response for zero degrees angle of attack is shown below in Figure 
3. The short period responses for the other angles of attack are not shown due to 


their similarity to the AOA=0 response. 
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Figure 3. Short Period Response for AOA = 0 degrees 


The initial conditions used to generate the short period response 
corresponded to a one-tenth scaling of the normalized complex eigenvector at 
=0. The normalized eigenvector was scaled in this manner to produce an initial 
perturbation quantity consistent with the magnitudes expected in flight. This 
technique took into account the phase relationship between each component of the 
eigenvector and produced the tue physical relationship between each component 
at any time. 

b. Root Locus 
It is interesting to note the changing position of the characteristic 


roots (eigenvalues) of the longitudinal system as angle of attack changes. The 
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| migration of the roots has a large impact on the character of the dynamic 
response in that the natural frequency and damping ratio are determined by the 
~~ Jocation of the roots in the complex plane. Shown below in Figure 4 is a root 


locus showing the migration of the short period roots. 

















l 


rad/sec 


rad/sec 





2 -1.5 =] -0.5 0 0.5 
Figure 4. Short Period Root Locus 


2. Dutch Roll Mode 

The F-14 displays an unstable dutch roll mode at high angle of attack 
due in part to degraded directional stability (Cag going positive to negative) 
which appears at around 15-20 degrees AOA, depending on aircraft speed and 
external stores loading. The Out of Control Flight Training Guide (Ref. 13] 
illustrates this degradation with speed and external loading as shown in Figures 5 
and 6. The variation of the directional stability parameter Cyg calculated from 
the data base is shown in Figure 7. The increase in dihedral effect (Cig getting 
more negative) and the decrease in roll damping (Cjp getting less negative) as 


. angle of attack increases also contribute significantly to the dutch roll instability. 
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Figure 5. Influence of Speed on Directional Stability 
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Figure 6. Influence of External Stores on Directional Stability 
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Figure 7. Variation of Ca8 with Angle of Attack from Database 


When this unstable mode is excited, the aircraft tends to diverge unless 
pilot inputs are made to control the motion. 
a. Time History Response 
The linear analysis reveals that the dutch roll mode is lightly 
damped at low angles of attack. This results in a slowly convergent oscillation 
for both zero and five degrees angle of attack. As angle of attack increases past 
ten degrees, the dutch roll mode goes unstable. At this point it is only slightly 
unstable, resulting in a very slow divergent oscillation. As angle of attack 
continues to increase, the dutch roll mode becomes extremely unstable, resulting 
in a rapid divergence. Figures 8, 9 and 10 show the dutch roll time history 


responses for AOA=0, 10 and 15 degrees, respectively. 
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Figure 8. Dutch Roll Response at AOA = 0 degrees 
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Figure 9. Dutch Roll Response at AOA = 10 degrees 
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Figure 10. Dutch Roll Response for AOA = 15 degrees 
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The initial conditions for all of the dutch roll responses used a one- 
tenth scaling of the normalized complex eigenvector, just as in the short period 
case. . 

b. Root locus 

An explanation for the behavior of the dutch roll mode is clear 
when an examination of the placement of the roots in the complex plane is made. 
A root locus is shown below in Figure 11 for the lateral-directional system's 
dutch roll roots. Notice the migration of the characteristic roots through the 
imaginary axis as the angle of attack approaches ten degrees. Once the roots 
migrate to the right hand half plane, the response changes from a convergent to a 


divergent oscillation. 


rad/sec 





4 
=0.4: 0.3 =0.2° 0.1 0 0.1 O02 O23 0.4 


Figure 11. Dutch Roll Root Locus 
As previously stated, the instability of the dutch roll mode at high 


angle of attack results primarily from changes in CnB, CIB and Clp. The root 


locus plots in Figures 12, 13 and 14 show the sensitivity of the roll, spiral and 
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dutch roll roots to individual changes in each of these parameters. The direction 
of root migration shown corresponds to variation of the stability parameter from 
~ the value at AOA=0 to AOA=25. All other parameters correspond to AOA=25 
values. When the changes in CnB, CIB and Clp occur simultaneously as AOA 
increases, the individual effects combine to produce a rapid onset of dutch roll 
“instability as shown previously in Figure 11. It was noted that had the values of 
Cn, CIB and Clp corresponding to zero AOA been held constant as AOA 


increased, the dutch roll roots would have been stable. 
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Figure 12. Root Locus for Cyg Variation 
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Figure 13. Root Locus for Cig Variation 
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Figure 14. Root Locus for Cjp Variation 
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B. NON-LINEAR ANALYSIS 

Now that an understanding of the effect of angle of attack on the linearized 
system dynamics has been gained, a comprehensive study of the nonlinearized 
_ system dynamics is in order. This begins by selecting the coupling option in the 
‘computer program, which has the effect of including the nonlinear terms in the 
equations of motion before numerical integration takes place. Additionally, 
initial conditions are chosen to demonstrate the effect of coupling from the 


| lateral-directional system to the longitudinal system. Specifically, the lateral- 


- directional initial conditions are the now familiar scaled eigenvector components, — 


while the longitudinal initial conditions are set to zero. 

At the low angles of attack where both linear modes are convergent, stable 
oscillations, one would expect that coupling of the two systems of equations 
would also produce a convergent, stable solution. This was exactly the case and 
is clearly evident from examination of Figure 15, the nonlinear response at zero 
degrees angle of attack. The changes in the lateral-directional parameters from 
the linear dutch roll response were imperceptible at this angle of attack. The 
effect of coupling on the longitudinal parameters is also evident from the plots, 
although the maximum amplitude of the longitudinal perturbations is on the 
order of a tenth of a degree. Although the longitudinal response shown here is 
imperceptible from a pilot's point of view, there is some significance to one 
feature which continues to appear at higher angles of attack. This is the 
frequency relationship between the roll angle response and the angle of attack 
response. In every case, the angle of attack perturbations occurred at twice the 


frequency of the roll response. This is somewhat intuitive, as angle of attack 
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Figure 15. Coupled Response at AOA = 0 degrees 
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perturbations can be expected to occur with either right or left wing movement. 


This tendency has been observed in the actual flight test time history traces 


~~ discussed earlier. 


A similar response was obtained for five degrees AOA. Once again, the 


nonlinear coupling produced very little change from the linear dutch roll 


-Tesponse and extremely small longitudinal perturbations. At ten degrees AOA, 


however, the nonlinear response failed to converge. Instead, a stable, constant 
magnitude oscillation developed in the lateral-directional parameters which is 
characteristic of a mild wing rock motion.. This response is shown in Figure 16. 
The maximum amplitudes of the longitudinal parameters due to the coupling are 
increasing, but are still insignificant in comparison to the overall response. Note 
that the mean value of the longitudinal perturbations in angle of attack and pitch 
rate are displaced from their respective equilibrium positions, while the pitch 
angle begins to fall off in an oscillatory manner. The presence of limit cycle 


behavior is clearly shown in this figure. 
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Figure 16. Coupled Response at AOA = 10 degrees 











Beyond ten degrees angle of attack, the coupling continued to produce 


oscillations in the lateral-directional parameters which were indicative of wing 


“*""*rock, while more dramatic changes began to appear in the longitudinal 


parameters. Responses for AOA=15, 20 and 25 degrees are shown on the 

following pages in Figures 17 through 19. 

_ The results of this portion of the analysis further demonstrated the 

development of wing rock limit cycles as a result of nonlinear coupling at high 

angle of attack. Although the maximum amplitudes for the roll angle and 

sideslip were somewhat higher than those given in NATOPS for the maneuver. 
flap configuration, the magnitudes were not unreasonable considering the 

approximations made to obtain the results. 

A few other important observations can be made by referring back to the 
plots of sideslip, roll angle, angle of attack and pitch angle for AOA=20, which 
have been enlarged in Figures 20a and 20bto show greater detail. The 
nonlinearity of the sideslip response near the maximum amplitude limits of the 
oscillation can be seen in Figure 20a. Instead of obtaining the smooth sinusoidal 
response normally associated with an underdamped system, this response 
illustrates that the coupling distorts the oscillation. Second, it is noted that as the 
perturbations in sideslip and roll angle begin, the frequency of the oscillations is 
very close to the dutch roll frequency. As the coupling takes effect and the limit 
cycle is established, the frequency of oscillation decreases by approximately 30 
percent. This trait was observed at all AOA's where a limit cycle was established 
and can be attributed to the nonlinear interaction of the longitudinal and lateral- 
directional modes. For the aircraft modeled in Ref. 4 which displayed a very 


mild dutch roll instability, this phenomenon was not ovserved. In that study, the 
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Figure 17. Coupled Response at AOA = 15 degrees 
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Figure 19. Coupled Response at AOA = 25 degrees 
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Figure 20b. Detailed View of Coupled Response at AOA = 20 degrees 
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frequencies of the natural modes were artificially controlled to maintain a 
harmonic relationship between the short period and the dutch roll response. 
‘Here, no such relationship exists. Further study may uncover the dependency of 
the limit cycle frequency on the characteristic frequencies of the natural modes. 
Finally, it was noted that as the limit cycle was established, the mean value of 
~ the angle of attack perturbation increased to a value which was slightly higher 
than its equilibrium position and the pitch angle began to drop. This tendency 
was also observed in the flight test results. Although the reasons for this are 
complex as a result of multiple dependency between the perturbation quantities in| 
the fully coupled nonlinear equations, one aspect of the wing rock motion may 
provide a clue as to its origin. Consider the lift equilibrium during the wing 
rock motion. The vertical component of lift at any time is proportional to the 
cosine of the roll angle, which changes continuously with time. Therefore, the 
average value of lift during one limit cycle oscillation is dependent upon the 
maximum amplitude of the wing rock and will always be lower than the 
equilibrium value of lift. As a result of the drop in the average lift per cycle, the 
aircraft starts to travel on a Curvelinear flight path. This response shows up as a 
loss in altitude. Additionally, the dependence of the nonlinear terms in the angle 
of attack and pitch angle equations ((15) and (17)) upon the average value of lift 
is readily apparent by the presence of the cos() term. This dependency may be 


the dominant feature which causes this type of response. 








VII. RESULTS 


The convergent nature of the linearized short period response and the 
divergent nature of the linearized dutch roll response for the same flight 
condition indicated that the numerical procedure used produced time history 
‘traces consistent with the behavior of stable and unstable linear systems, 
respectively. 

When the equations were modified to include the nonlinear terms, the 
responses for the low angles of attack (AOA=0 and 5 degrees) did not change 
appreciably, although the influence of the coupling was apparent in small 
longitudinal perturbations. At these angles of attack, both the longitudinal and 
lateral-directional characteristic roots were stable, resulting in an asymptotically 
stable solution in the sense of Liapunov [Ref. 14]. It has been demonstrated that 
although coupling between the two sets of equations occurred, the overall result 
failed to produce a limit cycle response. 

When the nonlinear equations were analyzed at the higher angles of attack, 
however, the presence of unstable dutch roll eigenvalues destroyed the 
asymptotic stability. Even so, the coupling between the longitudinal and lateral- 
directional equations still yielded a stable system response in the form of a wing 


rock limit cycle oscillation. 
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VIII. CONCLUSIONS 


A numerical an*lysis of the nenlinear equations of motion has been 
conducted to investigate the contributica of inertial coupling to the development 
“of wing rock in the F-14 aircraft. Actual wind tunnel data was used to develop 
all of the stability parameters for the analysis. Although a nurnber of 
simplifying assumptions were made, the analysis indicated that inertial coupling 
of a stable short period mode and an unstable dutch roll mode can result in a 
response very much like that encountered in the aircraft, especially for the - 
lateral-directional parameters. The trends displayed in the longitudinal 
parameters as a result of the coupling were consistent with flight tests; however, 
the magnitude of the excursions from the trim position for these parameters was 
greater than expected. Although these deviations appeared to be beyond reason, 
the results obtained here should not be discounted on this basis. The intention 
was to conduct a preliminary investigation into the mechanics of wing rock for 
the F-14 in hopes of uncovering a relatively simple explanation for the aircraft's 
behavior. In the pursuit of this goal, it has been demonstrated that a stable short 
period mode can feed damping energy into an unstable dutch roll mode to — 
produce a bounded wing rock type oscillation. Unquestionably there is a great 
deal of further research which must be done to fully understand the motion, 
including more detailed analyses which disregard some of the simplifications 
made in this analysis and which promise to yield results which are more 
consistent with the aircraft's actual response. The following section highlights a 


few of these techniques. 
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IX. RECOMMENDATIONS FOR FURTHER RESEARCH 


The analysis presented here is only a start. It serves to illustrate that 
numerical techniques can be used to solve the nonlinear equations of motion and 
to predict the response of aircraft subject to specified initial conditions. There 
exists a wide variety of related research topics that can further our understanding 
of nonlinear flight mechanics. Some of these topics are addressed in the 


paragraphs which follow. 


A. EIGHT DEGREE OF FREEDOM ANALYSIS 

Perhaps the least difficult of all sugges..ons to follow for continued research 
in this field would be the extension of the analysis to eight degrees of freedom. 
The eighth parameter of interest, perturbation velocity, was held constant during 
this analysis based on an approximation of constant velocity during the wing rock 
motion. To bring the analysis to eight degrees of freedom, the X-force equation 
would be introduced in its full non-linear form, the pertinent stability parameters 
obtained from further analysis of the available data base and the computer code 
modified. The influence of the additional degree of freedom on both the 
linearized analysis of the longitudinal parameters and the full non-linear analysis 
would warrant a full investigation of the changes in the dynamic response, 
especially given the tendency of the aircraft to end up in a nose down attitude 


where changes in velocity are sure to occur. 





B. TIME DEPENDENT STABILITY PARAMETER ANALYSIS 

A much more difficult undertaking would involve the incorporation of time 
dependent stability parameters into. the analysis. As stated previously, the data 
base contains an extensive list of stability parameters as a function of multiple 
variables; however, only derivative values corresponding to the steady, level 
flight trim condition were used for this investigation. Recall that as the aircraft 
moves, the pitch, roll and- yaw rates developed contribute to the stability of the 
aircraft and may play a significant role in the overall response of the aircraft. 
Furthermore, changes in AOA and sideslip angle as the aircraft moves about its 
equilibrium position also affect the stability parameters. Lastly, although this 
analysis was conducted controls fixed, .. > actual aircraft, control surface 
position also influences the stability c. «ic aircraft. Therefore, although this 
analysis proceeded under the approximation of constant stability characteristics 
as part of a planned buildup to a more detailed analysis, it nevertheless neglects 
the changing dynamic behavior of the aircraft as it translates and rotates. Ease 
of implementation for the first attempt at numerical creation of the wing rock 
motion from wind tunnel data was also a factor in choosing this approximation. 
Although this simplification may be appropriate for some aircraft operating at 
relatively low AOA, it may be the reason for obtaining such large perturbations 
in the longitudinal parameters at high AOA. 

I. order to implement such a change in the analysis, one of several options 
must be exercised to allow for variable stability parameters. One option might 
involve the development of subroutines which would allow the computer to 


conduct a table look up of dimensionless stability parameters and subsequent 
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conversion into dimensional form. This would require the availability of the 
“entire data base to the code, as well as restructuring of the code to accommodate 

~ ‘the table look up feature. re | : 

Another option might involve the use of multi-variable curve fitting to 

. obtain very accurate apEron manne of the ae stability parameters as 


- the aircraft moves. 


_C€. OPTIMIZATION OF NUMERICAL SOLUTION 

Still another possibility for future research in this area involves numerical 
optimization. Numerical techniques commonly used to obtain solutions for 
systems of differential equations vary widely in their accuracy and suitability for 
a given problem. The "exactness" of a numerical solution is very much 
dependent upon the methods used to numerically approximate the equations, the 
desired accuracy of the solution and the acceptable numerical cost in determining 
that solution. Similarly, the stability of a numerical procedure is also quite 
important in that it determines the acceptable variation of parameters, such as the 
time increment, which force a convergent solution. Care must be taken when 
conducting numerical analysis to account for the stability characteristics and 
desired accuracy of the solution. Therefore, a study attempting to optimize the 
numerical technique used to obtain the time history data points may be in order. 
For example, a number of different finite difference schemes could be used to 


obtain and compare solutions, computation time, overall efficiency, etc. 


D. INCORPORATION OF ACTUAL FLIGHT TEST RESULTS 
Once a numerical technique is developed to more accurately and efficiently 


model all aspects of the wing rock motion, the incorporation of additional flight 
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test results from an instrumented aircraft could be used to verify results obtained 
from the study. By utilizing the same initial conditions, aircraft configuration, 
“geometry and inertia characteristics as inputs into a similar analysis, it would be 


possible to try to numerically reconstruct the true response of the aircraft. 


‘E. NUMERICAL ANALYSIS OF F/A-18 WING ROCK 
‘In that a complete F/A-18 flight simulation program and data base are 
available for research use, a similar analysis conducted on that aircraft may 


provide some clues as to the mechanics behind F/A-18 wing rock. It is known 


that the aircraft displays a complex wing rock motion at very high AOA which 


cannot be damped out by pilot input, flight control computer or any combination 
of the two. The data base can be obtained from the NASA Ames Research 
Center/Dryden Flight Research Facility located at Edwards AFB, CA. 
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APPENDIX A.- DATA BASE MANIPULATION 


The stability derivatives needed to conduct this analysis were extracted from 


. tabular wind tunnel data provided by the Grumman Corporation (Ref. 9]. This 
_ appendix is provided to give the reader some insight as to how the tabular data 


"was used to extract the necessary information for the study. 


The first things needed to conduct the analysis at any particular angle of 


attack were the velocity, thrust required and control surface positions for steady, 


level flight (SLF) at that AOA. Assumptions of perfect lateral symmetry about 


the X-Z plane and the absence of gyroscopic effects led to the assignment of zero 
deflection for any lateral or directional control surfaces. The trim velocity and 
stabilator position were found by simultaneously solving the following two 


equations for alpha and 3s (stabilator position) in an iterative process: 
Sm eg = Sng * Sm thruse* Sm o% * ngs (23) 


CL= CLot Linrust * Lee + CL 5 5s (24) 


Specifically, the desired trim AOA was designated and a guess at the 
appropriate trim velocity was made. Cmcg is zero in SLF and Cy, can be 
calculated easily by using the assumed trim velocity and the weight of the 
aircraft. Walues for Cmo, Cma, Cio, and CLa were obtained by constructing 
crossplots of the data given in the tables to get the appropriate slope. Local 


derivatives were used in all cases where the data was nonlinear. Cmgs was taken 
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directly from the tables and C_g; was obtained by doing a coordinate 
transformation through the angle of attack of C,5, and Cz5, data in the table. In 
| ‘order to find Cyihruse and Cmihrusi, the ihiase: required for SLF had to be 
determined by finding the drag in SLF at the desired AOA. This was - 


accomplished by doing another coordinate transformation on the Cyo, Cro, Css 


and C,8, information in the table. A plot of Cp vs alpha was constructed from 

this data for each different stabilator setting given in the tables. From this plot, 

_an estimate of the drag cou!d be determined. This plot is shown below in Figure 
2 


Cp vs. Alpha 
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Figure 21. Cp vs Alpha 





With an estimate of the drag, thrust required is immediately known. This is 
converted into a pitching moment contribution and a lift contribution, which are 
. then converted into their final coefficient form. With all of the parameters for 
these two equations known, the equations are solved for alpha and $s. The initial 
_ assumed trim velocity is iterated until a match is made between the calculated 
AOA and the desired AOA. When a match is made, the stabilator setting 
calculated is used to update the determination of drag. The equations are solved 
again with this new drag figure and after only a few iterations, the calculation 
converges on the final result for stabilator setting. 

A convenient cross check for this procedure is available by building up a plot 
of Crcg vs. alpha by combining Cmo, Cmgs and Cmihrus: information. This 
procedure results in the plot shown below in Figure 22. This plot can be used to 
determine the stabilator setting required to maintain Cmcg = 0 (i.e., SLF) at any 
desired AOA. 

The determination of stabilator setting is important because some of the data 
in the tables is dependent on this parameter. Once this setting is known, all of 
the stability derivatives can be found. The following steps contain a brief 


descripticn of the remaining procedures to obtain these derivatives. 


1) Cyp, Cap, Cip and Czq are found by constructing crossplots of the data 
given in the tables and taking local derivatives. 
2) Cyp, Cyr, Cip, Cire Caps Car and Cg are taken directly from the tables. 


3) Ciq is found by doing a coordinate transformation of Cxg and Czq. 
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All of these dimensionless coefficients are then converted into dimensional 
form using standard conversion formulas. At this point, the derivatives are 


ready to be used in the analysis. 


Cmcg vs. Alpha 
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Figure 22. Cmcg vs Alpha 
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APPENDIX B.- COMPUTER PROGRAM 


100 REM F14 WING ROCK PROGRAM a, 
.-110 REM DATA FOR THE F14 AT AOA=20 DEGREES, Vtrim=213 fps aaa 
.120 REM A(, ) = PLANT MATRIX 

130 REM XK() = PERTURBATION QUANTITY (PQ) 

140 REM XK1() = ESTIMATION OF PQ AT NEXT FULL TIME STEP BASED 
150 REM ON DERIVATIVE AT CURRENT TIME STEP 

160 REM XK15( ) = ESTIMATION OF PQ AT HALF TIME STEP BASED ON 

170 REM DERIVATIVE AT CURRENT TIME STEP 

180 REM XK2() = ESTIMATION OF PQ AT FULL TIME STEP BASED ON 

190 REM DERIVATIVE AT HALF TIME STEP 

200 REM XSUM() = DERIVATIVE AT HALF OR FULL TIME STEP 

210 REM NL() = NONLINEAR DERIVATIVE TERMS 

220 REM XKMAX( ) = MAX VALUE OF PQ 


300 DIM A(10,10),XK(10),XK1(10),XK15(10),XK2(10) 
310 DIM XSUM(10),NL(10),XKMAX(10) 

320 ALPHA#=20 

330 VTRIM=213 

340 Pl=3.1415927# 

350 UsVTRIM*COS(ALPHA’PI/180) 

360 REM DEFINE INERTIAL PROPERTIES 

370 AlX=515091:AlY =232773!:AlZ=2756271:AlXZ=2654 
380 REM DEFINE TERMS IN EQNS OF MOTION 

390 AIQRL=(AIY-AIZ)/AIX ; AIPQN= (AIX-AIY)/AIZ 
400 AIRPMa(AIZ-AIX)/AlY 

410 AIPBM = -.16 : AZADOT#0 :‘THETO sALPHA*PI/180 
420 REM AIPBM is Malphadot, AZADOT=Zalphadot 

430 GU = 32.17/U 

440 DT=#.05 : DT2=.5°DT 


500 REM DATA FOR PLANTS. ROWS 1-4 LAT-DIR, ROWS 5-7 LONG 
510 DATA  -.0491,.0035,.1511,-1.0007,0,0,0 
520 DAT,. -6.338,-.5290,0,.6877,0,0,0 





_  §380 DATA 0,1,0,0,0,0,0 
540 DATA -.0515,-.0692,0,-.1186,0,0,0 
550 DATA 0,0,0,0,-.2671,.9659,-.0550 
“-- * * 3§80 DATA 0,0,0,0,-.2285,-.5278,.0088 
- 570 DATA 0,0,0,0,0,1,0 
' 580 REM READ IN PLANT DATA TO MATRIX A( ) 
590 FOR I=1 TO 7 
2600 FOR J=1TO7 
610 READ AiI,J) 
620 NEXT J 
630 NEXT |! 


700 REM INITIALIZE MATRICES FOR COMPUTATION 
710 FOR |=1 TO 7 

720 XK(I)=0 : XK1(I)=0 : XK15(I)eO : XK2(I)=0 
730 NEXT | 


800 REM DEFINE 1/10th SCALE DR EIGENVECTOR IC'S 

810 XK(1)= -.01434: XK(2)=0: XK(3)=-.09212 : XK(4)=.0055 
820 REM SET SP IC'S TO ZERO 

830 XK(5)= 0: XK(6)=0: XK(7)=0 


900 REM INPUT "USE INERTIAL COUPLING (Y/N)";A$ 
910 REM INPUT “CREATE DATA FILES (Y/N)":B$ 

920 A$="Y" : BS="Y" 

930 IF B$="Y" THEN GOSUB 2000 

940 REM PRINT INITIAL VALUES OF PQ'S 

950 GOSUB 2100 


1000 REM PERFORM NUMERICAL INTEGRATION 
1010 T=0 

1020 FOR K=1 TO 400 

1030 T=K*DT 


1100 REM CALCULATE DERIVATIVES AT CURRENT TIME STEP 
1100 REM USING PQ'S AT CURRENT TIME STEP 

1120 FOR Ilm1 TO 7 

1130 XSUM(I)=0 
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1140 
1150 
1160 


1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 


1280 


1290 
1300 
1310 
1320 


1400 
1410 
1420 
1430 
1440 
1450 
1460 


1500 
1510 
1520 
1530 
1540 
1550 
1560 
1570 
1580 
1590 
1600 


FOR J=1T07 
XSUM(I)=XSUM(I)+A(I JyPXK(d) 
NEXT J 


REM CALCULATE ESTIMATED PQ'S AT HALF AND FULL TIME 
REM STEP BASED ON DERIVATIVES AT CURRENT TIME STEP 
XK15(1)=XK(I)+DT2*XSUM(!) 
XK1(I)=XK(I)+DT*XSUM(I) 
NEXT | 
REM IF NECESSARY, CALCULATE NONLINEAR TERMS 
IF A$="N" THEN GOTO 1420 
GOSUB 2200 
REM CORRECT XK15() & XK1() FOR NONLINEAR TERMS 
FOR l-1 TO 7 
XK15(l)=XK15(1)+DT2*NL(I) 
XK1 (I)=XK1 (1)+DT*NL(I) 
NEXT | 


REM CALCUALTE DERIVATIVE AT HALF TIME STEP 
REM USING ESTIMATED PQ'S AT HALF TIME STEP 
FOR l=1 TO 7 
XSUM(I)=0 
FOR J=1TO7 
XSUM(I)=XSUM(I)+A(1,J)°XK15(J) 
NEXT J 


REM CALCULATE ESTIMATED PQ'S AT FULL TIME STEP 
REM BASED ON DERIVATIVE AT HALF TIME STEP 
XK2(N)=XK15(1)+OT2°XSUM(I) 
NEXT | 
REM IF NECESSARY, CALCULATE NONLINEAR TERMS 
IF A$="N" THEN GOTO 1720 
GOSUB 2400 
REM CORRECT XK2() FOR NONLINEAR TERMS 
FOR l= 1707 
XK2(I)=XK2(1)+DT2°NL(I) 
NEXT | 








1700 REM APPLY RICHARDSON'S EXTRAPOLATION TO GET 
1710 REM FINAL VALUE AT NEW TIME STEP 

1720 FOR le1 TO7 

“4730 = XK(I)=2*XK2(1)-XK1 (1) 


1800 REM FIND MAX VALUE OF EACH PQ 
1810 IF ABS(XK(I))>XKmax(il) THEN XKmax(i)-ABS(XK(I)) 


1900 NEXT I 
1910 REM PRINT RESULTS 
.1920 IF (INT(K/4)=(K/4)) THEN GOSUB 2100 
- 1930 NEXT K 
1940 IF B$="Y* THEN GOSUB 2700 
1990END 


2000 REM S/R FOR OPENING DATA FILES 
2010 OPEN "data" FOR OUTPUT AS #1 
2020 OPEN "“maxdata" FOR OUTPUT AS #2 
2030 RETURN 


2100 REM S/R PRINTS TIME AND STATE VECTOR ON SCREEN 
2110 REM AND STORES DATA IN FILE FOR SUBSEQUENT PLOTTING 
2120 PRINT USING “###.####°:T:XK(1);XK(2);XK(3):XK(4), 

| XK(5),XK(6),XK(7) 

2130 IF B$="Y" THEN GOSUB 2600 

2140 RETURN 


2200 REM S/R FOR FIRST USE OF NON LINEAR TERMS 
2210 THETA =THETO+XK(7) : CTHETA=COS(THETA) : 
TTHETA=TAN(THETA) 

2220 PHI=XK(3) : SPHI=SIN(PHI) : CPHI=COS(PHI) 

2230 NL(1)=XK(2)*XK(5)+GU*CTHETA*SPHI-A(1,3)*XK(3) 
2240 NL(2)=AIQRL*XK(6)*XK(4) 

2250 NL(3)=(XK(6)*SPHI+XK(4)*CPHI)*TTHETA 

2260 NL(4)=AIPQN*XK(2)*XK(6) 

2270 NL(5)=-XK(2)*XK(1)*U/(U-AZADOT)+32.17/(U-AZADOT) 
*(CTHETA*CPHI-COS(THETO)+SIN(THETO)*XK(7)) 
2280 NL(6)=AIRPM*XK(4)*XK(2)+AIPBM*NL(5) 
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2290 NL(7)=XK(6)*(CPHI- 1) XK(4)" SPHI 
2300 RETURN 


- 2400 REM S/R FOR SECOND USE OF NON LINEAR TERMS 
2410 THETA =THETO+XK15(7) : CTHETA=COS(THETA) : 
!  TTHETA»TAN(THETA) 
| 2420 PHleXK15(3) : SPHI=SIN(PHI) : CPHI=COS(PHI) . ~ 
ah 2490 NU(1)=XK15(2)°XK15(5)}+GU°CTHETA*SPHI- AC 3)°XK15(8)_ 
2 wee 4G0 ANL(2)mAIQRL*XK15(6)*XK15(4) “8 
2450 NL(3)=(XK15(6)* SPHI+XK15(4)*CPHI)* TTHETA ee 
2460 NL(4)=AIPQN*XK15(2)*XK15(6) ee ae 
2470 NL(5)=-XK15(2)*XK15(1)*U/(U-AZADOT)+32.17/(U-AZADOT) | 
*(CTHETA*CPHI-COS(THETO)+SIN(THETO)*XK(7)) 
2480 NL(6)=AIRPM*XK15(4)*XK15(2)+AIPBM*NL(5) 
72490 NL(7)=XK15(6)*(CPHI-1)-XK15(4)*SPHI 
2500 RETURN. 


t 
a 
it. 


_ .2600 REM S/R FOR STORING DATA TO FILE 
‘2610 PRINT #1, USING “###.# +8. #HHHH +H LHHHHE 
2620 REM CONTINUE PRINT USING +###.#### = +##H.HHHHH 
2630 REM CONTINUE +###.###0# +88. HHH +HHH AHHH, 
2640 REM CONTINUE T,XK(1),XK(2),XK(3),XK(4),XK(5),XK(6),XK(7) 
2650 RETURN 





2700 REM S/R TO CLOSE DATA FILE 

2710 CLOSE #1 

2720 PRINT #2, USING "+#.## +###.HHHHH = +HHH.HHHHH 
2730 REM CONTINUE PRINT USING +###.#####";zeta,XKmax(1) 
2740 REM CONTINUE ,XKmax(3),XKmax(5) 

2750 CLOSE #2 

2760 PRINT "DATA FILES CLOSED" ! 
2770 RETURN 
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